PIPELINE

For each sleep metric, we assess the GOF of the most complex model (including all predictors of interest) and we assess the following LMER models’ assumptions:

• independence between residuals and fitted values,

• the normality of both residuals’ and random effects’ distribution,

• that random effects are centered on zero,

TIB

Typically, TST shows a symmetrical shape that can be modeled with a normal distribution (as shown in the GRAPHICS section).

m <- lmer(TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)
# residuals normal fit (good)
plot(fitdist(resid(m),distr="norm",method="mle"))

# random effects normal fit (good)
plot(fitdist(ranef(m)$'ID'[[1]],distr="norm",method="mle"))

# LMER assumptions (good)
gridExtra::grid.arrange(plot(m,main="Resid vs Fitted"),
                        qqmath(residuals(m),main="Residual QQPlot"),
                        dotplot(ranef(m))[[1]],
                        plot(ranef(m),main="Random Effects QQplot")[[1]],
                        nrow=2)

Comments:

• The normal distribution shows a good fit on TIB’s residuals and random effects.

• All other LMER assumptions seem to be fulfilled.

Conclusions:

We can model TIB under the normality assumption.

TST

Like TIB, TST typically shows a symmetrical shape that can be modeled with a normal distribution (as shown in the GRAPHICS section).

m <- lmer(TST~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)
# residuals normal fit (good)
plot(fitdist(resid(m),distr="norm",method="mle"))

# random effects normal fit (good)
plot(fitdist(ranef(m)$'ID'[[1]],distr="norm",method="mle"))

# LMER assumptions (good)
gridExtra::grid.arrange(plot(m,main="Resid vs Fitted"),
                        qqmath(residuals(m),main="Residual QQPlot"),
                        dotplot(ranef(m))[[1]],
                        plot(ranef(m),main="Random Effects QQplot")[[1]],
                        nrow=2)

Comments:

• The normal distribution shows an excellent fit on TST’s residuals and a good fit on random effects.

• All other LMER assumptions seem to be fulfilled.

Conclusions:

We can model TST under the normality assumption.

SE

SE is bounded between 0 (when TST = 0) and 100 (when TST = TIB). Thus, a beta regression could be the best option, but at the time of these analyses no packages able to deal with mixed-effects beta regression were available (to our knowledge).

In our sample, only 4 on 82 showed an SE lower than 70%: EM07 (night 5), EM22 (night 36), EM36 (night 7) and EM53 (all 7 nights), whereas only one participant showed an SE = 99% and no participants showed a 100% SE.

# SE < 70%
sleep.long[!is.na(sleep.long$SE)&sleep.long$SE<70,
           c("ID","day","SE","TIB","TST","WASO","SOL")]

mean(sleep.long[sleep.long$ID=="MZ53","SE"])

## [1] 61.77

mean(sleep.long[sleep.long$ID=="MZ53","WASO"])

## [1] 148.2857

# SE = 99%
sleep.long[!is.na(sleep.long$SE)&sleep.long$SE==99,
           c("ID","day","SE","TIB","TST","WASO","SOL")]

The frequency distribution shows a quite normal shape, which is symmetric and centered approximately on 87-88%. The issue of modeling these data using a normal distribution is that the models would predict values above 100%, whereas TST cannot be higher than TIB (thus, SE cannot be higher than 100%).

hist(sleep.long$SE,breaks=100)

With no clues on how to model the data differently, we accept this limitation and we proceed by evaluating the fit of the normal distribution. Indeed, SE can be seen as a sort of transformation of TST, which can be reasonably modeled assuming residuals normality. SE results should confirm those obtained on TST, and we could focus our goal on assessing the impact of our predictors of interest on SE, instead of predicting meaningful SE values.

m <- lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)
# residuals normal fit (acceptable)
plot(fitdist(resid(m),distr="norm",method="mle"))

# random effects normal fit (acceptable)
plot(fitdist(ranef(m)$'ID'[[1]],distr="norm",method="mle"))

# LMER assumptions (acceptable)
gridExtra::grid.arrange(plot(m,main="Resid vs Fitted"),
                        qqmath(residuals(m),main="Residual QQPlot"),
                        dotplot(ranef(m))[[1]],
                        plot(ranef(m),main="Random Effects QQplot")[[1]],
                        nrow=2)

Comments:

• The normal distribution does not show a very good fit on SE’s residuals or random effects. This is especially due to the four outliers (EM07, EM22, EM36 and MZ53) showing SE lower than 70% (2+ hours of WASO, 2 to 5 hours of TST).

• All other LMER assumptions seem to be fulfilled (residuals are independent of fitted values, and random effects are normally distributed and centered on zero).

Let’s see what happens to residuals normality if we exclude the two outliers:

m <- lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),
          sleep.long[sleep.long$ID!="EM07"&sleep.long$ID!="EM22"&sleep.long$ID!="EM36"&sleep.long$ID!="MZ53",]) # excluding outliers
# residuals normal fit (acceptable)
plot(fitdist(resid(m),distr="norm",method="mle"))

Comments:

• the normal fit is much better without participants with SE <70%, and especially without EM07 and MZ53 (SE < 60%).

We can also try to log-transform the data considering the whole sample.

m <- lmer(log(SE)~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data=sleep.long)
# residuals normal fit (worse)
plot(fitdist(resid(m),distr="norm",method="mle"))

Comments:

• if we do not exclude the outliers, the log-transformation does not improve the fit of the normal distribution.

Conclusions:

Predicting SE scores assuming a normal residuals distribution implies that SE could be above 100% or (with much less probability) below 0%, which cannot be the case. However, with no clues on alternative ways to model our data, the normal distribution seems to be the better way. Looking at the data, we evaluate the normal fit on SE residuals as not excellent but at least acceptable. Removing four outliers would improve the fit.

SOL

SOL is a variable bounded on zero, typically skewed on the right, with few cases showing more than 1 hour of latency. Thus, we can try to log transform it to fit a model under normality, or we might use the Gamma family, which is bounded on zero.

# normal fit on log-transformed SOL (a constant is added to avoid cases of log(0))
m <- lmer(log(SOL+.05)~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)
# residuals normal fit (acceptable?)
plot(fitdist(resid(m),distr="norm",method="mle"))

# random effects normal fit (good)
plot(fitdist(ranef(m)$'ID'[[1]],distr="norm",method="mle"))

# LMER assumptions (acceptable)
gridExtra::grid.arrange(plot(m,main="Res. vs Fitted"),
                        qqmath(residuals(m),main="Residual QQPlot"),
                        dotplot(ranef(m))[[1]],
                        plot(ranef(m),main="Random Effects QQplot")[[1]],
                        nrow=2)

Comments:

• the fit of the normal distribution is not very good but acceptable

• all other LMER assumptions seem to be fulfilled (residuals are independent of fitted values, and random effects are normally distributed and centered on zero)

Then, we check the fit of the Gamma family.

m <- glmer(SOL+.05~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long,
           family="Gamma")
plot(fitdist(resid(m),distr="norm",method="mle"))

Conclusions:

The fit of the Gamma family on raw SOL data is worse than that of the normal distribution on log-transformed data. The best option seems to log-transform.

WASO

As TST, WASO usually shows a symmetrical shape that can be modeled with a normal distribution (as shown in the GRAPHICS section).

m <- lmer(WASO~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)
# residuals normal fit (good)
plot(fitdist(resid(m),distr="norm",method="mle"))

# random effects normal fit (good)
plot(fitdist(ranef(m)$'ID'[[1]],distr="norm",method="mle"))

# LMER assumptions (good)
gridExtra::grid.arrange(plot(m,main="Resid vs Fitted"),
                        qqmath(residuals(m),main="Residual QQPlot"),
                        dotplot(ranef(m))[[1]],
                        plot(ranef(m),main="Random Effects QQplot")[[1]],
                        nrow=2)

Comments:

• the normal distribution shows a good fit on WASO’s residuals and a good fit on random effects

• all other LMER assumptions seem to be fulfilled

Conclusions:

We can model WASO under the normality assumption.

BT

We expect to observe a symmetric shape also for BT.

m <- lmer(BT~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)
# residuals normal fit (acceptable)
plot(fitdist(resid(m),distr="norm",method="mle"))

# random effects normal fit (good)
plot(fitdist(ranef(m)$'ID'[[1]],distr="norm",method="mle"))

# LMER assumptions (good)
gridExtra::grid.arrange(plot(m,main="Resid vs Fitted"),
                        qqmath(residuals(m),main="Residual QQPlot"),
                        dotplot(ranef(m))[[1]],
                        plot(ranef(m),main="Random Effects QQplot")[[1]],
                        nrow=2)

# outliers
sleep.long[!is.na(sleep.long$BT)&sleep.long$BT>30,c("ID","day","BT","WT")]

Comments:

• the fit of the normal distribution is not very good but acceptable

• we can observe some outliers (especially EM22 and MZ36) showed BT higher than 34

• all other LMER assumptions seem to be fulfilled (residuals are independent of fitted values, and random effects are normally distributed and centered on zero)

Let’s see what happens if we exclude the two outliers:

m <- lmer(BT~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),
          sleep.long[sleep.long$ID!="EM22"&sleep.long$ID!="MZ36",]) # excluding outliers
# residuals normal fit (slighlty better)
plot(fitdist(resid(m),distr="norm",method="mle"))

Comments:

• the normal fit is slighlty better without participants EM22 and MZ36.

Also a log transformation slighlty improves the normal fit, but the improvement is not substantial.

m <- lmer(log(BT)~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)
# residuals normal fit (slighlty better)
plot(fitdist(resid(m),distr="norm",method="mle"))

Conclusions:

The normal fit on BT is not excellent but at least acceptable, removing two outliers seems to improve the fit. There is no evident reason to log transform the data.

WT

We expect to observe a symmetric shape also for WT.

m <- lmer(WT~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)
# residuals normal fit (good)
plot(fitdist(resid(m),distr="norm",method="mle"))

# random effects normal fit (good)
plot(fitdist(ranef(m)$'ID'[[1]],distr="norm",method="mle"))

# LMER assumptions (good)
gridExtra::grid.arrange(plot(m,main="Resid vs Fitted"),
                        qqmath(residuals(m),main="Residual QQPlot"),
                        dotplot(ranef(m))[[1]],
                        plot(ranef(m),main="Random Effects QQplot")[[1]],
                        nrow=2)

Comments:

• the normal distribution fits well both residuals and random effects for WT

• all other LMER assumptions seem to be fulfilled (residuals are independent of fitted values, and random effects are centered on zero)

Conclusions:

We can model WT under the normality assumption.

PIPELINE

For each sleep measure, we use the most complex model (with all predictors of interest) to proceed as follows:

• we compute the Cook’s distance by considering a cut-off set to 4/N,

• we re-estimate the parameters by removing all potentially influent cases identified with the CD, one-by-one,

• if a given case shows to substantially influence the parameters estimation, we replicate the influential analysis without that case,

• finally, we decide which cases can be removed.

TIB

# taking only ID and variables of interest (only continuous), specifying the model
data2check <- sleep.long[!is.na(sleep.long$TIB)&!is.na(sleep.long$week.time)&!is.na(sleep.long$SEX)&
                           !is.na(sleep.long$MEQ.r)&!is.na(sleep.long$PSQI)&!is.na(sleep.long$BDI.II)&
                           !is.na(sleep.long$ALCOL_week)&!is.na(sleep.long$FUMO..SI.NO.)&
                           !is.na(sleep.long$CAFFE_week),
                         c("ID","day","TIB","week.time","SEX","MEQ.r","PSQI","BDI.II","NAP",
                            "ALCOL_week","FUMO..SI.NO.","CAFFE_week")]
m <- lmer(TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data2check)

# computing Cook's dinstance and cutoff
cd <- cooks.distance(m)
data2check$cd <- cd
cutoff <- 4/length(levels(data2check$ID)) # 4/N

# selecting influential cases
index <- c(1:nrow(data2check)) # nrow
pm=cbind(data2check,cd,index)
pm=pm[order(pm$cd,decreasing = TRUE),] # order
pm$infl <- NA
for(i in 1:nrow(pm)){
  if(pm[i,"cd"] >= cutoff){ pm[i,"infl"] <- 1 } else { pm[i,"infl"] <- 0 }
}
PM <- pm

# Plotting Cook's distance
par(mfrow=c(1,1),mar=c(6, 5, 4, 2) + 0.1)
plot(cd, xlab="Index",ylab="Cook's distance",
     main=paste("Cook's Distance \ncut-off: 4/N =",round(cutoff,3),
                " ( NInfl =",nrow(pm[pm$infl==1,]),")"),
     pch=20,cex.lab=1.6,cex.axis=1.4,col.main="blue")
points(pm[pm$cd>=cutoff,"index"],
       pm[pm$cd>=cutoff,"cd"],pch=20,cex=1.3,col="red")
text(pm[pm$cd>=cutoff,"index"]-10,
     pm[pm$cd>=cutoff,"cd"]-0.005,
     paste(pm[pm$cd>=cutoff,"ID"],pm[pm$cd>=cutoff,"day"]),cex=.8)
abline(h=cutoff,
       lwd=1.8,col="blue",lty=2)

# participants showing cases above the cutoff
levels(as.factor(as.character(pm[pm$cd>cutoff,"ID"])))

##  [1] "EM07" "EM09" "EM22" "EM42" "em51" "MZ04" "MZ22" "MZ25" "MZ34" "MZ39"
## [11] "MZ40" "MZ46" "MZ55" "MZ56"

Comments:

• 21 observations (i.e., combination of participant and day) due to 14 participants are more extreme than others and show a CD above the cutoff.

• observations from MZ56 shows the highest Cook’s distance

In the next step, we re-estimate the parameters by excluding each potentially influential participant, one-by-one. We focus on substantial changes in the estimated parameters, such that the excluded cases implies a drastic change in the interpretation of the effects.

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(PM[PM$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- PM[PM$ID!=ID,]
  m4 <- lmer(TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

# week.time (the code for the remaining parameters is not showed)
p1 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,est))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")
p2 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,t))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")
gridExtra::grid.arrange(p1,p2,nrow=2)

Comments:

• Larger differences are shown by the parameters estimated for week.time (ranging from .20 to .28, with t-values ranging from 1.7 to 2.3) and Gender (ranging from -.39 to -.32, with t-values ranging from -.2.2 to -1.8).

• The most influential participants are EM42 (whose exclusion leads to a decrease of .04 in the parameter estimated for week.time, whose t-value decreases from 2.0 to 1.7), and MZ56 (whose exclusion leads to a decrease of .03 in the parameter estimated for Gender, whose t-value increases in absolute value from 1.95 to 2.2)

• None of the the other parameters seems to substantially change the interpretation of the results (i.e., non-significant effects remain non-significant, positive effects remain positive and negative effects remain negative).

Thus, we replicate the influential analysis by removing EM42 and MZ56 We start with the latter, which shows the highest Cook’s distance.

OUT <- data2check[data2check$ID!="MZ56",] # excluding influential case
m <- lmer(TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),OUT)

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(OUT[OUT$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- OUT[OUT$ID!=ID,]
  m4 <- lmer(TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

**Comments:*

• the exclusion of MZ56 increases the effect of Gender, which becomes more consistent across the sample.

• the role of week.time is still inconsistent, with participant EM42 showing the strongest influence (leading to an increase of .04 in the estimated parameter, with t-value increasing from 1.5 to 1.9).

• none of the other parameters shows to be substantially influenced by the exclusion of MZ56.

Thus, we replicate the influential analysis by removing EM42

OUT <- OUT[OUT$ID!="EM42",] # excluding influential case
m <- lmer(TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),OUT)

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(OUT[OUT$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- OUT[OUT$ID!=ID,]
  m4 <- lmer(TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

Comments:

• the exclusion of EM42 decreases the effect of Week time, which becomes more consistent across the sample.

• none of the other parameters shows to be substantially influenced by the exclusion of EM42 (i.e., non-significant effects remain non-significant, positive effects remain positive and negative effects remain negative).

Conclusions:

We decide to perform the analyses on TIB without MZ56 and EM42.

TST

# taking only ID and variables of interest (only continuous), specifying the model
data2check <- sleep.long[!is.na(sleep.long$TST)&!is.na(sleep.long$week.time)&!is.na(sleep.long$SEX)&
                           !is.na(sleep.long$MEQ.r)&!is.na(sleep.long$PSQI)&!is.na(sleep.long$BDI.II)&
                           !is.na(sleep.long$ALCOL_week)&!is.na(sleep.long$FUMO..SI.NO.)&
                           !is.na(sleep.long$CAFFE_week),
                         c("ID","day","TST","week.time","SEX","MEQ.r","PSQI","BDI.II","NAP",
                            "ALCOL_week","FUMO..SI.NO.","CAFFE_week")]
m <- lmer(TST~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data2check)

# computing Cook's dinstance and cutoff
cd <- cooks.distance(m)
data2check$cd <- cd
cutoff <- 4/length(levels(data2check$ID)) # 4/N

# selecting influential cases
index <- c(1:nrow(data2check)) # nrow
pm=cbind(data2check,cd,index)
pm=pm[order(pm$cd,decreasing = TRUE),] # order
pm$infl <- NA
for(i in 1:nrow(pm)){
  if(pm[i,"cd"] >= cutoff){ pm[i,"infl"] <- 1 } else { pm[i,"infl"] <- 0 }
}
PM <- pm

# Plotting Cook's distance
par(mfrow=c(1,1),mar=c(6, 5, 4, 2) + 0.1)
plot(cd, xlab="Index",ylab="Cook's distance",
     main=paste("Cook's Distance \ncut-off: 4/N =",round(cutoff,3),
                " ( NInfl =",nrow(pm[pm$infl==1,]),")"),
     pch=20,cex.lab=1.6,cex.axis=1.4,col.main="blue")
points(pm[pm$cd>=cutoff,"index"],
       pm[pm$cd>=cutoff,"cd"],pch=20,cex=1.3,col="red")
text(pm[pm$cd>=cutoff,"index"]-10,
     pm[pm$cd>=cutoff,"cd"]-0.005,
     paste(pm[pm$cd>=cutoff,"ID"],pm[pm$cd>=cutoff,"day"]),cex=.8)
abline(h=cutoff,
       lwd=1.8,col="blue",lty=2)

# participants showing cases above the cutoff
levels(as.factor(as.character(pm[pm$cd>cutoff,"ID"])))

##  [1] "EM06" "EM07" "EM09" "EM22" "EM23" "EM36" "EM42" "em51" "EM53" "MZ04"
## [11] "MZ08" "MZ22" "MZ25" "MZ34" "MZ39" "MZ40" "MZ44" "MZ46" "MZ50" "MZ55"
## [21] "MZ56"

Comments:

• 30 observations (i.e., combination of participant and day) due to 21 participants are more extreme than others and show a CD above the cutoff

In the next step, we re-estimate the parameters by excluding each potentially influential participant, one-by-one. We focus on substantial changes in the estimated parameters, such that the excluded cases implies a drastic change in the interpretation of the effects.

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(PM[PM$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- PM[PM$ID!=ID,]
  m4 <- lmer(TST~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

# week.time (the code for the remaining parameters is not showed)
p1 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,est))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")
p2 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,t))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")
gridExtra::grid.arrange(p1,p2,nrow=2)

Comments:

• Larger differences are shown by the parameters estimated for week.time (ranging from 12 to 16, with t-values ranging from 2.0 to 2.4), MEQ-R (ranging from -1.25 to 0.00, with t-values ranging from -1.00 to 0.00, and MZ44 and MZ56 leading to an average decrease of .4 in the estimated parameter), ALCOHOL (ranging from -1.6 to -.2, with t-values ranging from -1.0 to 0.0, and MZ44 leading to an absolute increase of .75 in the estimated parameter), SMOKE (ranging from -15 to -8, with t-values ranging from -1.4 to -.8, and EM06, EM22 and EM36 leading to an average absolute increase of .2 in the estimated parameter).

• EM22 and EM36 are among the most influential cases. Interestingly, they are also two of the four outliers found for SE.

• None of the potentially influential cases seems to substantially change the interpretation of the results (i.e., non-significant effects remain non-significant, positive effects remain positive and negative effects remain negative).

Conclusions:

We decide to keep all participants in TST analyses.

SE

# taking only ID and variables of interest (only continuous), specifying the model
data2check <- sleep.long[!is.na(sleep.long$SE)&!is.na(sleep.long$week.time)&!is.na(sleep.long$SEX)&
                           !is.na(sleep.long$MEQ.r)&!is.na(sleep.long$PSQI)&!is.na(sleep.long$BDI.II)&
                           !is.na(sleep.long$ALCOL_week)&!is.na(sleep.long$FUMO..SI.NO.)&
                           !is.na(sleep.long$CAFFE_week),
                         c("ID","day","SE","week.time","SEX","MEQ.r","PSQI","BDI.II","NAP",
                            "ALCOL_week","FUMO..SI.NO.","CAFFE_week")]
m <- lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data2check)

# computing Cook's dinstance and cutoff
cd <- cooks.distance(m)
data2check$cd <- cd
cutoff <- 4/length(levels(data2check$ID)) # 4/N

# selecting influential cases
index <- c(1:nrow(data2check)) # nrow
pm=cbind(data2check,cd,index)
pm=pm[order(pm$cd,decreasing = TRUE),] # order
pm$infl <- NA
for(i in 1:nrow(pm)){
  if(pm[i,"cd"] >= cutoff){ pm[i,"infl"] <- 1 } else { pm[i,"infl"] <- 0 }
}
PM <- pm

# Plotting mahalanobis distance
par(mfrow=c(1,1),mar=c(6, 5, 4, 2) + 0.1)
plot(cd, xlab="Index",ylab="Cook's distance",
     main=paste("Cook's Distance \ncut-off: 4/N =",round(cutoff,3),
                " ( NInfl =",nrow(pm[pm$infl==1,]),")"),
     pch=20,cex.lab=1.6,cex.axis=1.4,col.main="blue")
points(pm[pm$cd>=cutoff,"index"],
       pm[pm$cd>=cutoff,"cd"],pch=20,cex=1.3,col="red")
text(pm[pm$cd>=cutoff,"index"]-10,
     pm[pm$cd>=cutoff,"cd"]-0.005,
     paste(pm[pm$cd>=cutoff,"ID"],pm[pm$cd>=cutoff,"day"]),cex=.8)
abline(h=cutoff,
       lwd=1.8,col="blue",lty=2)

# participants showing cases above the cutoff
levels(as.factor(as.character(pm[pm$cd>cutoff,"ID"])))

##  [1] "EM02" "EM04" "EM07" "EM08" "EM10" "EM12" "EM16" "EM22" "EM23" "EM29"
## [11] "EM36" "EM39" "EM54" "MZ03" "MZ04" "MZ05" "MZ09" "MZ15" "MZ21" "MZ25"
## [21] "MZ30" "MZ35" "MZ36" "MZ38" "MZ53" "MZ56"

Comments:

• 39 observations (i.e., combination of participant and day) due to 26 participants are more extreme than others and show a CD above the cutoff

In the next step, we re-estimate the parameters by excluding each potentially influential participant, one-by-one. We focus on substantial changes in the estimated parameters, such that the excluded cases implies a drastic change in the interpretation of the effects.

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(data2check[data2check$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- data2check[data2check$ID!=ID,]
  m4 <- lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

# week.time (the code for the remaining parameters is not showed)
p1 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,est))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  geom_text(data=parameters[parameters$pred=="week.timewe" & (parameters$est<=.1|parameters$est>=.3),],
            aes(label=ID),nudge_x=-1) + theme(axis.text.x = element_blank())
p2 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,t))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  geom_text(data=parameters[parameters$pred=="week.timewe" & (parameters$est<=.1|parameters$est>=.3),],
            aes(label=ID),nudge_x=-1) + theme(axis.text.x = element_blank())
gridExtra::grid.arrange(p1,p2,nrow=2)

Comments:

• Larger differences are shown by the parameter estimated for NAP (ranging from -.6 to -1.1, with t-values ranging from -1.6 to -2.3, and EM07 showing the strongest influence, leading to an increase of .4 in the estimated NAP coefficient). Further differences are shown by the parameters estimated for sex (ranging from -2.0 to -1.0, with t-values ranging from -1.6 to -1.0, and MZ53 leading to an absolute increase of .75 in the estimated parameter), and smoke (ranging from -3 to -1.5, with t-values ranging from -2.5 to -1.5, and MZ53 leading to an absolute increase of 1.0 in the estimated parameter).

• MZ53 shows a degree of influence on almost all parameters, although the interpretation of the results does not change substantially (i.e., non-significant effects remain non-significant, positive effects remain positive and negative effects remain negative). The only exception is smoke, whose effect seems to strongly depend on MZ53.

• MZ53 is also one of the four outliers found for SE, with all nights showing an SE < .70.

Thus, we replicate the influential analysis by removing EM07 and MZ53. We start with the former.

OUT <- data2check[data2check$ID!="EM07",] # excluding influential case
m <- lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),OUT)

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(OUT[OUT$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- OUT[OUT$ID!=ID,]
  m4 <- lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

Comments:

• The exclusion of EM07 decreases the effect of NAP and makes the estimated parameter more consistent.

• Participant MZ53 continues to show a degree of influence on many parameters, and particularly on smoke.

Thus, we replicate the influential analysis by removing MZ53.

OUT <- OUT[OUT$ID!="MZ53",] # excluding influential case
m <- lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),OUT)

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(OUT[OUT$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- OUT[OUT$ID!=ID,]
  m4 <- lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

Comments:

• The exclusion of MZ53 decreases the effect of smoke and makes the estimated parameter more consistent. None of the other parameters show to be substantially affected by the removal of this participant.

• We do not observe other potentially influential cases.

Conclusions:

We decide to perform the analyses on SE without EM07 and MZ53.

SOL

# taking only ID and variables of interest (only continuous), specifying the model
data2check <- sleep.long[!is.na(sleep.long$SOL)&!is.na(sleep.long$week.time)&!is.na(sleep.long$SEX)&
                           !is.na(sleep.long$MEQ.r)&!is.na(sleep.long$PSQI)&!is.na(sleep.long$BDI.II)&
                           !is.na(sleep.long$ALCOL_week)&!is.na(sleep.long$FUMO..SI.NO.)&
                           !is.na(sleep.long$CAFFE_week),
                         c("ID","day","SOL","week.time","SEX","MEQ.r","PSQI","BDI.II","NAP",
                            "ALCOL_week","FUMO..SI.NO.","CAFFE_week")]
m <- lmer(log(SOL+.05)~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),sleep.long)

# computing Cook's dinstance and cutoff
cd <- cooks.distance(m)
data2check$cd <- cd
cutoff <- 4/length(levels(data2check$ID)) # 4/N

# selecting influential cases
index <- c(1:nrow(data2check)) # nrow
pm=cbind(data2check,cd,index)
pm=pm[order(pm$cd,decreasing = TRUE),] # order
pm$infl <- NA
for(i in 1:nrow(pm)){
  if(pm[i,"cd"] >= cutoff){ pm[i,"infl"] <- 1 } else { pm[i,"infl"] <- 0 }
}
PM <- pm

# Plotting mahalanobis distance
par(mfrow=c(1,1),mar=c(6, 5, 4, 2) + 0.1)
plot(cd, xlab="Index",ylab="Cook's distance",
     main=paste("Cook's Distance \ncut-off: 4/N =",round(cutoff,3),
                " ( NInfl =",nrow(pm[pm$infl==1,]),")"),
     pch=20,cex.lab=1.6,cex.axis=1.4,col.main="blue")
points(pm[pm$cd>=cutoff,"index"],
       pm[pm$cd>=cutoff,"cd"],pch=20,cex=1.3,col="red")
text(pm[pm$cd>=cutoff,"index"]-10,
     pm[pm$cd>=cutoff,"cd"]-0.005,
     paste(pm[pm$cd>=cutoff,"ID"],pm[pm$cd>=cutoff,"day"]),cex=.8)
abline(h=cutoff,
       lwd=1.8,col="blue",lty=2)

# participants showing cases above the cutoff
levels(as.factor(as.character(pm[pm$cd>cutoff,"ID"])))

##  [1] "EM04" "EM05" "EM08" "EM09" "EM12" "EM18" "EM23" "EM24" "EM26" "EM44"
## [11] "EM53" "EM54" "MZ02" "MZ04" "MZ05" "MZ14" "MZ15" "MZ21" "MZ22" "MZ23"
## [21] "MZ31" "MZ35" "MZ36" "MZ39" "MZ47" "MZ49" "MZ50" "MZ51" "MZ55" "MZ57"

Comments:

• 38 observations (i.e., combination of participant and day) due to 30 participants are more extreme than others and show a CD above the cutoff

In the next step, we re-estimate the parameters by excluding each potentially influential participant, one-by-one. We focus on substantial changes in the estimated parameters, such that the excluded cases implies a drastic change in the interpretation of the effects.

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(data2check[data2check$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- data2check[data2check$ID!=ID,]
  m4 <- lmer(log(SOL+.05)~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)

  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

# week.time (the code for the remaining parameters is not showed)
p1 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,est))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  geom_text(data=parameters[parameters$pred=="week.timewe" & parameters$est<=(-.1),],
            aes(label=ID),nudge_x=-1) + theme(axis.text.x = element_blank())
p2 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,t))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  geom_text(data=parameters[parameters$pred=="week.timewe" & parameters$est<=(-.1),],
            aes(label=ID),nudge_x=-1) + theme(axis.text.x = element_blank())
gridExtra::grid.arrange(p1,p2,nrow=2)

Comments:

• None of the estimated parameters seems to be substantially affected by potentially influential cases.

• Two participants show stronger influence on several parameters: MZ15 for sex, MEQ-r, PSQI, and BDI; MZ23 for alcohol and coffee. However, also these participants do not change the interpretation of the results (i.e., non-significant effects remain non-significant, positive effects remain positive and negative effects remain negative).

Conclusions:

We decide to keep all participants in SOL analyses.

WASO

# taking only ID and variables of interest (only continuous), specifying the model
data2check <- sleep.long[!is.na(sleep.long$WASO)&!is.na(sleep.long$week.time)&!is.na(sleep.long$SEX)&
                           !is.na(sleep.long$MEQ.r)&!is.na(sleep.long$PSQI)&!is.na(sleep.long$BDI.II)&
                           !is.na(sleep.long$ALCOL_week)&!is.na(sleep.long$FUMO..SI.NO.)&
                           !is.na(sleep.long$CAFFE_week),
                         c("ID","day","WASO","week.time","SEX","MEQ.r","PSQI","BDI.II","NAP",
                            "ALCOL_week","FUMO..SI.NO.","CAFFE_week")]
m <- lmer(WASO~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data2check)

# computing Cook's dinstance and cutoff
cd <- cooks.distance(m)
data2check$cd <- cd
cutoff <- 4/length(levels(data2check$ID)) # 4/N

# selecting influential cases
index <- c(1:nrow(data2check)) # nrow
pm=cbind(data2check,cd,index)
pm=pm[order(pm$cd,decreasing = TRUE),] # order
pm$infl <- NA
for(i in 1:nrow(pm)){
  if(pm[i,"cd"] >= cutoff){ pm[i,"infl"] <- 1 } else { pm[i,"infl"] <- 0 }
}
PM <- pm

# Plotting mahalanobis distance
par(mfrow=c(1,1),mar=c(6, 5, 4, 2) + 0.1)
plot(cd, xlab="Index",ylab="Cook's distance",
     main=paste("Cook's Distance \ncut-off: 4/N =",round(cutoff,3),
                " ( NInfl =",nrow(pm[pm$infl==1,]),")"),
     pch=20,cex.lab=1.6,cex.axis=1.4,col.main="blue")
points(pm[pm$cd>=cutoff,"index"],
       pm[pm$cd>=cutoff,"cd"],pch=20,cex=1.3,col="red")
text(pm[pm$cd>=cutoff,"index"]-10,
     pm[pm$cd>=cutoff,"cd"]-0.005,
     paste(pm[pm$cd>=cutoff,"ID"],pm[pm$cd>=cutoff,"day"]),cex=.8)
abline(h=cutoff,
       lwd=1.8,col="blue",lty=2)

# participants showing cases above the cutoff
levels(as.factor(as.character(pm[pm$cd>cutoff,"ID"])))

##  [1] "EM04" "EM07" "EM09" "EM10" "EM12" "EM17" "EM22" "EM29" "EM38" "EM39"
## [11] "EM50" "EM53" "EM54" "EM55" "MZ03" "MZ16" "MZ17" "MZ25" "MZ30" "MZ34"
## [21] "MZ36" "MZ42" "MZ46" "MZ49" "MZ53" "MZ56"

Comments:

• 44 observations (i.e., combination of participant and day) due to 26 participants are more extreme than others and show a CD above the cutoff

In the next step, we re-estimate the parameters by excluding each potentially influential participant, one-by-one. We focus on substantial changes in the estimated parameters, such that the excluded cases implies a drastic change in the interpretation of the effects.

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(PM[PM$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- PM[PM$ID!=ID,]
  m4 <- lmer(WASO~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

# week.time (the code for the remaining parameters is not showed)
p1 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,est))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  geom_text(data=parameters[parameters$pred=="week.timewe" & parameters$est>=1.5,],
            aes(label=ID),nudge_x=-1) + theme(axis.text.x = element_blank())
p2 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,t))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  geom_text(data=parameters[parameters$pred=="week.timewe" & parameters$est>=1.5,],
            aes(label=ID),nudge_x=-1) + theme(axis.text.x = element_blank())
gridExtra::grid.arrange(p1,p2,nrow=2)

Comments:

• Larger differences are shown by the parameter estimated for smoke (ranging from 7 to 13, with t-values ranging from 1.4 to 2.2), whose effect is strongly influenced by MZ53 (leading to an increase of 4 in the estimated smoke coefficient).

• MZ53 is the participant influencing more coefficients (week.time, sex, MEQ-r, PSQI, alcohol, and smoke), followed by EM54 (affecting the parameter for alcohol), MZ16 (influencing the effect of coffee) and MZ25 (influencing the effect of BDI-II). MZ53 is the most influential case (this participant is also one of the outliers for SE).

Thus, we replicate the influential analysis by removing MZ53.

OUT <- data2check[data2check$ID!="MZ53",] # excluding influential case
m <- lmer(WASO~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),OUT)

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(OUT[OUT$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- OUT[OUT$ID!=ID,]
  m4 <- lmer(WASO~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

Comments:

• The exclusion of MZ53 decreases and makes more consistent the parameter estimated for smoke. Slight changes are observable also for other parameters.

• There are still some participants showing a degree of influence on the estimation of the parameters, especially MZ25 (affecting sex and PSQI), MZ25 (affecting PSQI, BDI-II, and smoke), and MZ36 (affecting sex, alcohol, and smoke). However, the exclusion of these participant does not substantially change the estimated parameters.

Conclusions:

We decide to exclude only **MZ53* from WASo analyses.

BT

# taking only ID and variables of interest (only continuous), specifying the model
data2check <- sleep.long[!is.na(sleep.long$BT)&!is.na(sleep.long$week.time)&!is.na(sleep.long$SEX)&
                           !is.na(sleep.long$MEQ.r)&!is.na(sleep.long$PSQI)&!is.na(sleep.long$BDI.II)&
                           !is.na(sleep.long$ALCOL_week)&!is.na(sleep.long$FUMO..SI.NO.)&
                           !is.na(sleep.long$CAFFE_week),
                         c("ID","day","BT","week.time","SEX","MEQ.r","PSQI","BDI.II","NAP",
                            "ALCOL_week","FUMO..SI.NO.","CAFFE_week")]
m <- lmer(BT~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data2check)

# computing Cook's dinstance and cutoff
cd <- cooks.distance(m)
data2check$cd <- cd
cutoff <- 4/length(levels(data2check$ID)) # 4/N

# selecting influential cases
index <- c(1:nrow(data2check)) # nrow
pm=cbind(data2check,cd,index)
pm=pm[order(pm$cd,decreasing = TRUE),] # order
pm$infl <- NA
for(i in 1:nrow(pm)){
  if(pm[i,"cd"] >= cutoff){ pm[i,"infl"] <- 1 } else { pm[i,"infl"] <- 0 }
}
PM <- pm

# Plotting mahalanobis distance
par(mfrow=c(1,1),mar=c(6, 5, 4, 2) + 0.1)
plot(cd, xlab="Index",ylab="Cook's distance",
     main=paste("Cook's Distance \ncut-off: 4/N =",round(cutoff,3),
                " ( NInfl =",nrow(pm[pm$infl==1,]),")"),
     pch=20,cex.lab=1.6,cex.axis=1.4,col.main="blue")
points(pm[pm$cd>=cutoff,"index"],
       pm[pm$cd>=cutoff,"cd"],pch=20,cex=1.3,col="red")
text(pm[pm$cd>=cutoff,"index"]-10,
     pm[pm$cd>=cutoff,"cd"]-0.005,
     paste(pm[pm$cd>=cutoff,"ID"],pm[pm$cd>=cutoff,"day"]),cex=.8)
abline(h=cutoff,
       lwd=1.8,col="blue",lty=2)

# participants showing cases above the cutoff
levels(as.factor(as.character(pm[pm$cd>cutoff,"ID"])))

##  [1] "EM07" "EM08" "EM09" "EM22" "EM23" "EM29" "EM36" "EM40" "EM42" "EM45"
## [11] "MZ03" "MZ04" "MZ17" "MZ22" "MZ25" "MZ34" "MZ40" "MZ50" "MZ53" "MZ54"
## [21] "MZ55" "MZ56"

Comments:

• 28 observations (i.e., combination of participant and day) due to 22 participants are more extreme than others and show a CD above the cutoff

In the next step, we re-estimate the parameters by excluding each potentially influential participant, one-by-one. We focus on substantial changes in the estimated parameters, such that the excluded cases implies a drastic change in the interpretation of the effects.

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(data2check[data2check$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- PM[PM$ID!=ID,]
  m4 <- lmer(BT~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

# week.time (the code for the remaining parameters is not showed)
p1 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,est))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  theme(axis.text.x = element_blank())
p2 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,t))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  theme(axis.text.x = element_blank())
gridExtra::grid.arrange(p1,p2,nrow=2)

Comments:

• Two participants show stronger influence on several parameters: EM36 for week.time, sex, MEQ-r, PSQI, and smoke; MZ22 for EQr, alcohol and coffee. Alcohol is the mostly influenced parameter, whose magitude is decreased .072 (t = 1.89) to .047 (t = 1.22).

• However, the exclusion of these participants does not change the interpretation of the results (i.e., non-significant effects remain non-significant, positive effects remain positive and negative effects remain negative).

Conclusions:

We decide to keep all participants in BT data analysis.

WT

# taking only ID and variables of interest (only continuous), specifying the model
data2check <- sleep.long[!is.na(sleep.long$WT.dayAfter)&!is.na(sleep.long$week.time)&!is.na(sleep.long$SEX)&
                           !is.na(sleep.long$MEQ.r)&!is.na(sleep.long$PSQI)&!is.na(sleep.long$BDI.II)&
                           !is.na(sleep.long$ALCOL_week)&!is.na(sleep.long$FUMO..SI.NO.)&
                           !is.na(sleep.long$CAFFE_week),
                         c("ID","day","WT.dayAfter","week.time","SEX","MEQ.r","PSQI","BDI.II","NAP",
                            "ALCOL_week","FUMO..SI.NO.","CAFFE_week")]
m <- lmer(WT.dayAfter~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data2check)

# computing Cook's dinstance and cutoff
cd <- cooks.distance(m)
data2check$cd <- cd
cutoff <- 4/length(levels(data2check$ID)) # 4/N

# selecting influential cases
index <- c(1:nrow(data2check)) # nrow
pm=cbind(data2check,cd,index)
pm=pm[order(pm$cd,decreasing = TRUE),] # order
pm$infl <- NA
for(i in 1:nrow(pm)){
  if(pm[i,"cd"] >= cutoff){ pm[i,"infl"] <- 1 } else { pm[i,"infl"] <- 0 }
}
PM <- pm

# Plotting mahalanobis distance
par(mfrow=c(1,1),mar=c(6, 5, 4, 2) + 0.1)
plot(cd, xlab="Index",ylab="Cook's distance",
     main=paste("Cook's Distance \ncut-off: 4/N =",round(cutoff,3),
                " ( NInfl =",nrow(pm[pm$infl==1,]),")"),
     pch=20,cex.lab=1.6,cex.axis=1.4,col.main="blue")
points(pm[pm$cd>=cutoff,"index"],
       pm[pm$cd>=cutoff,"cd"],pch=20,cex=1.3,col="red")
text(pm[pm$cd>=cutoff,"index"]-10,
     pm[pm$cd>=cutoff,"cd"]-0.005,
     paste(pm[pm$cd>=cutoff,"ID"],pm[pm$cd>=cutoff,"day"]),cex=.8)
abline(h=cutoff,
       lwd=1.8,col="blue",lty=2)

# participants showing cases above the cutoff
levels(as.factor(as.character(pm[pm$cd>cutoff,"ID"])))

##  [1] "EM01" "EM07" "EM08" "EM09" "EM17" "EM22" "EM30" "EM36" "EM39" "EM47"
## [11] "em51" "MZ03" "MZ04" "MZ08" "MZ17" "MZ22" "MZ25" "MZ39" "MZ40" "MZ49"
## [21] "MZ54" "MZ56"

Comments:

• 32 observations (i.e., combination of participant and day) due to 22 participants are more extreme than others and show a CD above the cutoff

In the next step, we re-estimate the parameters by excluding each potentially influential participant, one-by-one. We focus on substantial changes in the estimated parameters, such that the excluded cases implies a drastic change in the interpretation of the effects.

# parameters with all participants
parameters <- as.data.frame(round(summary(m)$coefficients,3))
parameters$ID <- "All"
IDs <- levels(as.factor(as.character(data2check[data2check$cd>cutoff,"ID"]))) # ID-day couples with CD > cutoff
for(ID in IDs){
  new.data <- PM[PM$ID!=ID,]
  m4 <- lmer(WT.dayAfter~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),new.data)
  new.data <- as.data.frame(round(summary(m4)$coefficients,3))
  new.data$ID <- ID
  parameters <- rbind(parameters,new.data)
}
parameters <- cbind(parameters,rep(c("Intercept","week.timewe","SEXM","MEQ.r","PSQI","BDI.II","NAP",
                                     "ALCOL_week","FUMO..SI.NO.1","CAFFE_week"),
                                   length(IDs)+1))
colnames(parameters) <- c("est","st.err","t","ID","pred")
parameters

# week.time (the code for the remaining parameters is not showed)
p1 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,est))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  theme(axis.text.x = element_blank())
p2 <- ggplot(parameters[parameters$pred=="week.timewe",],aes(ID,t))+
  geom_point()+ggtitle("week.time effect - original (red) vs. removing cases")+xlab("")+
  geom_point(data=parameters[parameters$pred=="week.timewe"&parameters$ID=="All",],colour="red")+
  theme(axis.text.x = element_blank())
gridExtra::grid.arrange(p1,p2,nrow=2)

Comments:

• None of the estimated parameters seems to be substantially affected by potentially influential cases.

• MZ22 is the participant influencing more parameters, and particoularly MEQ-r, alcohol, smoke, and coffee comsumption. However, the exclusion of these participants does not change the interpretation of the results (i.e., non-significant effects remain non-significant, positive effects remain positive and negative effects remain negative). The strongest influence of this participant is on alcohol assumption, whose estimated parameter changes from .05 (t = 1.82) to .03 (t = 1.08) after the removal of MZ22.

infl <- influence.ME::influence(m,"ID")
cd <- influence.ME::cooks.distance.estex(infl, parameter=8, sort=TRUE)[,1]
cd[which(names(cd)=="MZ22")]

##      MZ22 
## 0.5399679

4/nlevels(data2check$ID) # cutoff

## [1] 0.04878049

Conclusions:

We decide to keep all participants in WT data analysis.

PIPELINE

For each sleep measure:

• We specify the models including the predictors in the order mentioned above (see Aims and Content).

• The set of models is compared using the LRT and the AIC. The AIC weight (AICw) is used to express AIC in a probabilistic way (i.e., each model is associated with a probability, from 0 to 1, to be the most plausible model).

• The parameters estimated by the selected model are reported with their 95% confidence intervals (CI) obtained using bootstrap with the percentile method (with 10,000 replicates), their standard error (SE), and t-value.

The function modelComp is used to optimize the process and generate the output.

modComp <- function(data=sleep.long,measure="TIB",influentials=c("MZ56","EM42"),
                    CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2){
  
  require(lme4); require(MuMIn)
  
  # excluding influential cases
  df <- data
  for(i in 1:length(influentials)){ df <- df[df$ID!=influentials[i],] }
  
  # model specification
  models <- list(lmer(formula=gsub("TIB",measure,
                                   "TIB~(1|ID)"), df), # null model
                 lmer(formula=gsub("TIB",measure,
                                   "TIB~week.time+(1|ID)"),df), # week time
                 lmer(formula=gsub("TIB",measure,
                                   "TIB~week.time+SEX+(1|ID)"), df), # SEX
                 lmer(formula=gsub("TIB",measure,
                                   "TIB~week.time+SEX+MEQ.r+(1|ID)"), df), # MEQ.r
                 lmer(formula=gsub("TIB",measure,
                                   "TIB~week.time+SEX+MEQ.r+PSQI+(1|ID)"), df), # PSQI
                 lmer(formula=gsub("TIB",measure,
                                   "TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+(1|ID)"), df), # BDI
                 lmer(formula=gsub("TIB",measure,
                                   "TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+(1|ID)"), df), # NAP
                 lmer(formula=gsub("TIB",measure,
                                   "TIB~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID)"), df)) # sub
  
  # likelihood ratio test
  lrt <- as.data.frame(anova(models[[1]],models[[2]],models[[3]],models[[4]],models[[5]],
                             models[[6]],models[[7]],models[[8]]))
  mc <- cbind(model=c("null","week.time","SEX","MEQ.r","PSQI","BDI.II","NAP","sub"),
              round(lrt[,c(1,2,5:7)],digits),p=round(lrt$`Pr(>Chisq)`,digits+1))
  
  # AIC weight
  mc$AICw <- round(Weights(mc$AIC),digits)
  
  # model selection (based on LRT and AIC)
  bestMod <- mc[which(mc$AICw==max(mc$AICw)),"model"]
  cat("selected model:",as.character(bestMod),"\n\n")
  print(summary(models[[which(mc$model==bestMod)]]))
  
  # selected model's parameters
  coeff <- round(summary(models[[which(mc$model==bestMod)]])$coefficients,digits)
  
  # boostrap CI
  coeff <- cbind(coeff,round(confint.merMod(models[[which(mc$model==bestMod)]],
                                            parm=3:(length(fixef(models[[which(mc$model==bestMod)]]))+2),
                                            level=CI.level,method=CI.method,boot.type=boot.type,nsim=nsim),digits))
  
   # joining model comparison and coefficients info
  if(nrow(coeff) < nrow(mc)){
    empty.row <- rep("",nrow(mc)-nrow(coeff))
    out <- cbind(mc[,c("model","Chisq","Chi Df","p","AICw")],
                 rbind(coeff[,c(1:2,4:5,3)],
                       cbind(Estimate=empty.row,'Std. Error'=empty.row,
                             '2.5 %'=empty.row,'97.5 %'=empty.row,'t value'=empty.row)))
    } else if(nrow(coeff) > nrow(mc)){
      empty.row <- rep("",nrow(coeff)-nrow(mc))
      out <- cbind(rbind(mc[,c("model","Chisq","Chi Df","p","AICw")],
                         cbind(model=empty.row,Chisq=empty.row,'Chi Df'=empty.row,
                               p=empty.row,AICw=empty.row)),
                   coeff[,c(1:2,4:5,3)]) }
  rownames(out) <- 1:nrow(out)
  out[,c("Estimate","Std. Error","2.5 %","97.5 %","t value")] <- lapply(out[,c("Estimate","Std. Error","2.5 %","97.5 %","t value")],
                                                                        as.character)
  out[,c("Estimate","Std. Error","2.5 %","97.5 %","t value")] <- lapply(out[,c("Estimate","Std. Error","2.5 %","97.5 %","t value")],
                                                                        as.numeric)
  
  return(out)
}

TIB

TIB is modeled assuming a normal distribution for residuals and excluding participant MZ56 and EM42 (see previous sections).

(TIB <- modComp(data=sleep.long,measure="TIB",influentials=c("MZ56","EM42"),
                CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2))

## selected model: SEX 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: TIB ~ week.time + SEX + (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 1885.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1865 -0.5921  0.0757  0.5594  3.6386 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.2214   0.4706  
##  Residual             1.5526   1.2460  
## Number of obs: 557, groups:  ID, 80
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   7.7926     0.1080  72.151
## week.timewe   0.1775     0.1167   1.521
## SEXM         -0.4755     0.1493  -3.184
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw
## week.timewe -0.309       
## SEXM        -0.654 -0.003

Comments:

• The model including the effect of both Week time and Gender is the only one showing a significantly higher logLikelihood compared to the null model for TIB.

• The model including Week time and Gender is also the model with the lowest AIC, followed by the model including MEQ.r (AICw = .21).

• TIB is only predicted by participants’ Gender (with males showing shorter sleep time compared to females).

Below, the effects are plotted and a summary table is generated.

# effect plot
plot(effects::allEffects(lmer(TIB~week.time+SEX+(1|ID),sleep.long)))

Let’s see what the results would have been if MZ56 and EM42 were not excluded.

modComp(data=sleep.long,measure="TIB",influentials="none",
        CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2)

## selected model: SEX 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: TIB ~ week.time + SEX + (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 1980.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.4029 -0.5651  0.0680  0.5462  3.5227 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.2317   0.4813  
##  Residual             1.6921   1.3008  
## Number of obs: 571, groups:  ID, 82
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   7.7155     0.1093  70.606
## week.timewe   0.2392     0.1203   1.988
## SEXM         -0.4162     0.1526  -2.727
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw
## week.timewe -0.315       
## SEXM        -0.644 -0.003

Comments:

• Including MZ56 and EM42 would have led to an higher parameter estimated for Week Time, with higher TIB during weekend compared to weekdays (coeff. = .24, SE = .12, t = 1.99), while mantaining the effect of Gender.

• Also in this case, the model with both Week time and Gender would have been selected based on AIC.

Conclusions:

TIB is significantly predicted only by participants’ gender, with shorter TIB in males compared to females. The effect of Week time was less evident and due to participant EM42.

TST

TST is modeled assuming a normal distribution for residuals and including all participants (see previous sections).

(TST <- modComp(data=sleep.long,measure="TST",influentials="none",
                CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2))

## selected model: SEX 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: TST ~ week.time + SEX + (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 1850.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.3662 -0.5428  0.0680  0.5632  3.8629 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.2774   0.5267  
##  Residual             1.3010   1.1406  
## Number of obs: 571, groups:  ID, 82
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   6.7188     0.1069  62.822
## week.timewe   0.2305     0.1055   2.184
## SEXM         -0.4875     0.1509  -3.229
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw
## week.timewe -0.282       
## SEXM        -0.651 -0.003

Comments:

• The model including the effect of week.time, and that including also sex are the only models associated to a significantly higher logLikelihood compared to the null model for TST.

• The model including only sex and week.time is also the model with the lowest AIC, follo. The dAIC suggests that the former is about 2.5 times more plausible than the latter.

• TST is predicted by participants’ sex (with males showing shorter sleep duration), and by week time (with weekend being associated to longer sleep duration than weekdays).

Below, the effects are plotted.

# effect plot
plot(effects::allEffects(lmer(TST~week.time+SEX+(1|ID),sleep.long)))

Conclusions:

Total sleep time is only predicted by participants’ gender (with males showing a shorter duration than females) and week time (with longer TST during weekend compared to week days).

SE

SE is modeled assuming a normal distribution for residuals and excluding participant EM07 and MZ53 (see previous sections).

(SE <- modComp(data=sleep.long,measure="SE",influentials=c("EM07","MZ53"),
               CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2))

## selected model: null 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: SE ~ (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 3045.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8514 -0.4831  0.0708  0.5593  3.3593 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 13.596   3.687   
##  Residual              9.912   3.148   
## Number of obs: 557, groups:  ID, 80
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  86.8377     0.4333   200.4

Comments:

• None of the specified models shows a significantly higher logLikelihood than the null model.

• Consistently, the null model shows the lowest AIC (AICw = .37), with model 1 (including only week.time) showing the second lowest AIC (AICw = .25).

Let’s see what the results would have been if EM07 and MZ53 were not excluded.

(SE <- modComp(data=sleep.long,measure="SE",influentials="none",
               CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2))

## selected model: null 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: SE ~ (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 3247.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.1059 -0.4471  0.0718  0.5300  3.1171 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 21.56    4.644   
##  Residual             11.91    3.451   
## Number of obs: 571, groups:  ID, 82
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  86.4205     0.5328   162.2

summary(lmer(SE~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data=sleep.long))

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## SE ~ week.time + SEX + MEQ.r + PSQI + BDI.II + NAP + ALCOL_week +  
##     FUMO..SI.NO. + CAFFE_week + (1 | ID)
##    Data: sleep.long
## 
## REML criterion at convergence: 3238.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -7.9580 -0.4446  0.0905  0.5389  3.0974 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 20.07    4.480   
##  Residual             11.84    3.441   
## Number of obs: 571, groups:  ID, 82
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)   91.040487   2.621406  34.730
## week.timewe    0.165767   0.318890   0.520
## SEXM          -1.786462   1.170683  -1.526
## MEQ.r         -0.129762   0.151364  -0.857
## PSQI          -0.003708   0.231190  -0.016
## BDI.II        -0.044624   0.082999  -0.538
## NAP1          -0.911895   0.454523  -2.006
## ALCOL_week     0.122842   0.177051   0.694
## FUMO..SI.NO.1 -2.611608   1.244778  -2.098
## CAFFE_week    -0.117683   0.069559  -1.692
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw SEXM   MEQ.r  PSQI   BDI.II NAP1   ALCOL_ FUMO..
## week.timewe -0.036                                                        
## SEXM        -0.212 -0.001                                                 
## MEQ.r       -0.809  0.000  0.028                                          
## PSQI        -0.403 -0.001  0.151  0.039                                   
## BDI.II      -0.124  0.000 -0.149  0.140 -0.547                            
## NAP1        -0.019  0.056  0.019 -0.001  0.001 -0.016                     
## ALCOL_week  -0.031  0.000 -0.445  0.026 -0.130  0.185 -0.010              
## FUMO..SI.NO -0.168 -0.001  0.044  0.163  0.062 -0.095 -0.018 -0.291       
## CAFFE_week  -0.042 -0.001  0.197 -0.212  0.028 -0.054 -0.016 -0.332 -0.134

Comments:

• Including EM07 and MZ53 would have led the LRT to indicate the most complex model as the best one.

• In contrast, the null model is still the model with the highest AICw.

• The most complex model shows an effect of NAP and Smoke, with lower SE on days with NAP compared to days without NAP (coef. = -.91, SE = .45, t = -2.01) and on smokers compared to non-smokers (coef. = -2.61, SE = 1.18, t = -2.10).

Conclusions:

None of the included predictors exerted a significant effect on SE. The effect of NAP was due to participant EM07, and the effect of Smoke was due to participant MZ53.

SOL

SOL is log-transformed and modeled assuming a normal distribution for residuals. All participants are included in the analysis (see previous sections).

(SOL <- modComp(data=sleep.long,measure="log(SOL+.05)",influentials="none",
                CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2))

## selected model: sub 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## log(SOL + 0.05) ~ week.time + SEX + MEQ.r + PSQI + BDI.II + NAP +  
##     ALCOL_week + FUMO..SI.NO. + CAFFE_week + (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 2228.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0754 -0.3678  0.1508  0.6414  2.1997 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.6508   0.8067  
##  Residual             2.3802   1.5428  
## Number of obs: 571, groups:  ID, 82
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    1.39856    0.56185   2.489
## week.timewe   -0.04922    0.14293  -0.344
## SEXM           0.60455    0.25030   2.415
## MEQ.r         -0.04750    0.03233  -1.469
## PSQI          -0.03430    0.04937  -0.695
## BDI.II         0.03627    0.01774   2.045
## NAP1           0.23561    0.19593   1.203
## ALCOL_week    -0.08856    0.03780  -2.343
## FUMO..SI.NO.1 -0.13342    0.26593  -0.502
## CAFFE_week     0.03043    0.01485   2.049
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw SEXM   MEQ.r  PSQI   BDI.II NAP1   ALCOL_ FUMO..
## week.timewe -0.075                                                        
## SEXM        -0.211 -0.001                                                 
## MEQ.r       -0.807  0.000  0.027                                          
## PSQI        -0.403 -0.001  0.151  0.040                                   
## BDI.II      -0.124  0.001 -0.152  0.140 -0.546                            
## NAP1        -0.038  0.054  0.039 -0.001  0.003 -0.032                     
## ALCOL_week  -0.031  0.001 -0.446  0.026 -0.130  0.187 -0.020              
## FUMO..SI.NO -0.168 -0.002  0.043  0.164  0.063 -0.094 -0.037 -0.290       
## CAFFE_week  -0.041 -0.002  0.196 -0.212  0.028 -0.053 -0.032 -0.332 -0.133

Comments:

• The model including the effect of BDI-II and the most complex model (with substances) are the only models associated with a significantly higher logLikelihood compared to the null model for SE.

• Coherently, the last model shows the lowest AIC (AICw = .44), with the model with BDI-II showing the second lowest AIC (AICw = .17).

• SOL is mainly predicted by participants’ sex (with longer SOL for males), BDI-II score, and the weekly amount of consumed alcohol and coffee.

Below, the effects are plotted.

# effect plots
plot(effects::allEffects(lmer(log(SOL + 0.05) ~ week.time + SEX + MEQ.r + PSQI + BDI.II + NAP +  
                                ALCOL_week + FUMO..SI.NO. + CAFFE_week + (1 | ID),sleep.long))[c(2,5,7,9)])

Conclusions:

The model comparison suggested a positive relationship between SOL and depressive symptoms, a negative relationship between SOL and the weekly amount of alcoholic units, a positive relationship between SOL and the weekly amount of coffees, and a higher SOL among smokers compared to non-smokers. Moreover, males showed shorted SOL than females.

WASO

WASO is modeled assuming a normal distribution for residuals and by excluding participant MZ53 (see previous sections).

(WASO <- modComp(data=sleep.long,measure="WASO",influentials="MZ53",
                 CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2))

## selected model: null 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: WASO ~ (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 5040.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4436 -0.5891 -0.0494  0.5612  4.3741 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 285.1    16.88   
##  Residual             340.9    18.46   
## Number of obs: 564, groups:  ID, 81
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   51.710      2.031   25.46

Comments:

• None of the specified models shows a significantly higher logLikelihood than the null model.

• Consistently, the null model shows the lowest AIC (AICw = .46), with model 1 (including only week.time) showing the second lowest AIC (AICw = .32).

Let’s see what would have been the results if MZ53 was not excluded.

modComp(data=sleep.long,measure="WASO",influentials="none",
        CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2)

## selected model: null 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: WASO ~ (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 5146.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4217 -0.5571 -0.0399  0.5367  4.3105 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 392.8    19.82   
##  Residual             355.1    18.85   
## Number of obs: 571, groups:  ID, 82
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   52.885      2.327   22.73

Comments:

• although in section 6 we observed and increased effect of Smoke due to participant MZ53, the inclusion of this participant does not change the results of the model comparison.

Conclusions:

None of the included predictors exerted a significant effect on WASO. The effect of Smoke was not substantial, and mainly due to participant MZ53.

BT

BT is modeled assuming a normal distribution for residuals and including all participants (see previous sections).

(BT <- modComp(data=sleep.long,measure="BT",influentials="none",
               CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2))

## selected model: MEQ.r 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: BT ~ week.time + SEX + MEQ.r + (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 2062
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3704 -0.5531 -0.0790  0.4277  7.4683 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.8255   0.9086  
##  Residual             1.7353   1.3173  
## Number of obs: 571, groups:  ID, 82
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 26.77508    0.45815  58.442
## week.timewe  0.68769    0.12187   5.643
## SEXM         0.42605    0.23026   1.850
## MEQ.r       -0.14951    0.03185  -4.694
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw SEXM  
## week.timewe -0.076              
## SEXM        -0.300 -0.002       
## MEQ.r       -0.937  0.000  0.073

# summary of the most complex model
summary(lmer(BT~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data=sleep.long))

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## BT ~ week.time + SEX + MEQ.r + PSQI + BDI.II + NAP + ALCOL_week +  
##     FUMO..SI.NO. + CAFFE_week + (1 | ID)
##    Data: sleep.long
## 
## REML criterion at convergence: 2076.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.3333 -0.5553 -0.0620  0.4410  7.4006 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.7522   0.8673  
##  Residual             1.7383   1.3185  
## Number of obs: 571, groups:  ID, 82
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)   26.358165   0.563701  46.759
## week.timewe    0.688531   0.122158   5.636
## SEXM           0.182421   0.251344   0.726
## MEQ.r         -0.143271   0.032475  -4.412
## PSQI           0.002539   0.049594   0.051
## BDI.II         0.001892   0.017814   0.106
## NAP1           0.010741   0.169745   0.063
## ALCOL_week     0.071955   0.037978   1.895
## FUMO..SI.NO.1  0.421426   0.267109   1.578
## CAFFE_week     0.002939   0.014918   0.197
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw SEXM   MEQ.r  PSQI   BDI.II NAP1   ALCOL_ FUMO..
## week.timewe -0.064                                                        
## SEXM        -0.211 -0.001                                                 
## MEQ.r       -0.808  0.000  0.027                                          
## PSQI        -0.403 -0.001  0.151  0.040                                   
## BDI.II      -0.124  0.001 -0.151  0.140 -0.547                            
## NAP1        -0.033  0.055  0.034 -0.001  0.002 -0.028                     
## ALCOL_week  -0.031  0.001 -0.445  0.026 -0.130  0.186 -0.017              
## FUMO..SI.NO -0.168 -0.002  0.043  0.164  0.063 -0.094 -0.032 -0.290       
## CAFFE_week  -0.041 -0.002  0.196 -0.212  0.028 -0.054 -0.028 -0.332 -0.133

Comments:

• The model including the effect of week.time shows a significantly higher logLikelihood compared to the null model for BT. Moreover, the LRT suggests a significant contribution of MEQ-r score over week time and sex, and a significant contribution of the inclusion of substances over the other predictors.

• In contrast, the model with MEQ-r is selected as the most plausible model based on the AIC (AICw = .44), with the most complex model showing the second best AIC (AICw = .31).

• In the selected model, BT is predicted by the week time, with weekend days showing later BT compared to week days. Moreover, BT was predicted by the MEQ-r score, with earlier BT being predicted by higher scores (i.e., indicating morning circadian preferences). The effect of sex is not substantial.

• The most complex model (selected based on the LRT but not the AIC), does not show any substantial effect of substances.

Below, the effects are plotted.

plot(effects::allEffects(lmer(BT ~ week.time + SEX + MEQ.r + (1|ID),sleep.long))[c(1,3)])

Conclusions:

The model’s comparison suggested a negative relationship between BT and the score on the MEQ-r, and later BT during the weekend compared to weekdays.

WT

WT is modeled assuming a normal distribution for residuals and including all participants (see previous sections). Note that for testing the effects of week.time on this variable, we should use WT.dayAfter, which is referred to the subsequent day (e.g., WT on day 7 is Monday WT). In this way, Monday to Friday are considered weekdays, whereas Saturday and Sunday are considered weekend days.

(WT <- modComp(data=sleep.long,measure="WT.dayAfter",influentials="none",
               CI.method="boot",boot.type="perc",CI.level=.95,nsim=10000,digits=2))

## selected model: MEQ.r 
## 
## Linear mixed model fit by REML ['lmerMod']
## Formula: WT.dayAfter ~ week.time + SEX + MEQ.r + (1 | ID)
##    Data: df
## 
## REML criterion at convergence: 1891.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.9263 -0.6002 -0.0830  0.5683  4.2610 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.4798   0.6927  
##  Residual             1.3376   1.1566  
## Number of obs: 569, groups:  ID, 82
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 10.462207   0.362748  28.842
## week.timewe  0.890752   0.107298   8.302
## SEXM        -0.002575   0.182162  -0.014
## MEQ.r       -0.146897   0.025204  -5.828
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw SEXM  
## week.timewe -0.081              
## SEXM        -0.300 -0.003       
## MEQ.r       -0.936 -0.003  0.074

# summary of the most complex model
summary(lmer(WT~week.time+SEX+MEQ.r+PSQI+BDI.II+NAP+ALCOL_week+FUMO..SI.NO.+CAFFE_week+(1|ID),data=sleep.long))

## Linear mixed model fit by REML ['lmerMod']
## Formula: 
## WT ~ week.time + SEX + MEQ.r + PSQI + BDI.II + NAP + ALCOL_week +  
##     FUMO..SI.NO. + CAFFE_week + (1 | ID)
##    Data: sleep.long
## 
## REML criterion at convergence: 1967
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8816 -0.5948 -0.0227  0.5381  3.7976 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID       (Intercept) 0.4098   0.6402  
##  Residual             1.5125   1.2298  
## Number of obs: 569, groups:  ID, 82
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)   10.430786   0.446613  23.355
## week.timewe    0.252725   0.114017   2.217
## SEXM          -0.178386   0.199107  -0.896
## MEQ.r         -0.143042   0.025708  -5.564
## PSQI           0.005645   0.039253   0.144
## BDI.II        -0.008849   0.014102  -0.628
## NAP1           0.099003   0.156954   0.631
## ALCOL_week     0.053742   0.030052   1.788
## FUMO..SI.NO.1  0.371894   0.211621   1.757
## CAFFE_week    -0.001581   0.011803  -0.134
## 
## Correlation of Fixed Effects:
##             (Intr) wk.tmw SEXM   MEQ.r  PSQI   BDI.II NAP1   ALCOL_ FUMO..
## week.timewe -0.075                                                        
## SEXM        -0.210  0.001                                                 
## MEQ.r       -0.807  0.000  0.027                                          
## PSQI        -0.403 -0.002  0.150  0.040                                   
## BDI.II      -0.124  0.001 -0.152  0.140 -0.545                            
## NAP1        -0.038  0.052  0.037 -0.004  0.005 -0.030                     
## ALCOL_week  -0.031  0.000 -0.446  0.026 -0.129  0.187 -0.018              
## FUMO..SI.NO -0.169 -0.004  0.041  0.164  0.064 -0.093 -0.036 -0.289       
## CAFFE_week  -0.041 -0.002  0.196 -0.211  0.027 -0.054 -0.033 -0.332 -0.133

Comments:

• The model including the effect of week.time shows a significantly higher logLikelihood compared to the null model for BT. Moreover, the LRT suggests a significant contribution of the MEQ-r score over week time and sex, and a significant contribution of substances’ effect over the other predictors.

• In contrast, the model with MEQ-r is selected as the most plausible model based on the AIC (AICw = .40), with the most complex model showing the second best AIC (AICw = .35).

• In the selected model, WT is predicted by the week time, with weekend days showing later BT compared to week days. Moreover, BT is predicted by the MEQ-r score, with earlier BT being predicted by higher scores (i.e., indicating morning circadian preferences).

• The most complex model (selected based on the LRT but not the AIC), does not show any substantial effect of substances.

Below, the effects are plotted.

plot(effects::allEffects(lmer(WT ~ week.time + SEX + MEQ.r + (1|ID),sleep.long))[c(1,3)])

Conclusions:

The model comparison suggested a negative relationship between WT and the score on the MEQ-r, and later WT during the weekend compared to weekdays.

FINAL TABLE

Here, the outputs of the model comparisons are combined to generate the results table.

(results <- rbind(TIB,TST,SE,SOL,WASO,BT,WT))

write.csv(results,"ModelComp.csv")